Perturbation theory as a new analytical approach to the MEG forward problem 
for realistic volume conductor modeling of the head 

G. Nolte 1 , T. Fieseler 2 , and G. Curio 1 

1 Neurophysics Group, Dept, of Neurology, Klinikum Benjamin Franklin, Berlin, Germany; 

2 Institute of Medicine, Research Center Jiilich, Germany. 


1 Introduction 

Analytical solutions to the magnetic forward prob¬ 
lem are known only for special volume conductors 
[1, 2, 3], which are insufficient to describe detailed 
head shapes. The forward calculation for a realistic 
volume conductor can so far only be solved by numer¬ 
ical algorithms, such as the boundary element method 
(BEM) [4, 5, 6, 7, 8, 9], which pose severe problems 
for superficial sources due to the discretization of the 
volume conductor surface. Here, we present an ana¬ 
lytical approach using perturbation theory where we 
exploit the fact that for realistic head shapes the devia¬ 
tion of the inner skull from a spherical approximation 
is small compared to its radius. 


2 Theory 


2.1 Description of surface 


In order to find an approximate analytical solution we 
have to analytically define the surface [10, 11]. We 
assume that a function / exists such that the surface is 
given by the image of a function G : [0. ir] x [0. 2ir] —> 
iR 3 with 
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R is the unperturbed, constant radius of the spherical 
approximation of the real volume conductor. Here we 
expand / in the basis of spherical harmonics up to 
order P 
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The (f3 P q) are now the parameters which define the 
surface. Though the expansion can in principle be 
made arbitrarily accurate large spatial frequencies re¬ 
quire larger computational effort (see below): a trun¬ 
cation of the series at not too large values for P 



Figure 1: Surfaces 


should lead to a sufficient approximation of the vol¬ 
ume surface. Here we chose P = 6 corresponding 
to (P + l) 2 = 49 parameters. We found that e.g. 
the surface corresponding to the liquor-skull bound¬ 
ary can be parametrized in this way with an average 
error of about 1.5 mm. In Fig. 1 we show the fit result 
in comparison to the original surface. Note, that this 
approach leads to a perfectly smooth surface descrip¬ 
tion. 


2.2 Formal forward calculation 

Expressing the external magnetic field as a surface in¬ 
tegral leads to a linear relationship between magnetic 
field and surface potential. At first sight it seems that 
hence an n.th order perturbative solution of the mag¬ 
netic field requires the solution of the corresponding 
electric problem [10]. However, calculating the total 
magnetic field from the radial component (with re¬ 
spect to the unperturbed sphere) leads to a shift of 
the orders because a secondary current can only con¬ 
tribute to the radial component of B if also its direc¬ 
tion is non-radial and is hence absorbing at least one 
order of the perturbation. 

The starting point of the theoretical analysis is the 




standard integral relating the external magnetic field 
to an integral over the surface dC 
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where V is the surface potential. 

The measure dS can be made explicit by standard 
methods of vector analysis leading to 

dS = d@d$r sin© (reR - , 
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Insertion of (5) into (4) leads to a complicated expres¬ 
sion for B. Interestingly, we found the following sim¬ 
ple and yet exact integral for the contribution of the 
secondary currents to the radial component of-B: 


are fixed numbers which can be calculated in closed 
form for all indices. 

Splitting the whole computation in an initialization 
and an actual forward calculation for each dipole is 
recommended if a large number of sources is consid¬ 
ered. The initialization now consists of two steps 

1. For the surface coefficients /3p q which parame¬ 
terizes the perturbation around the sphere with 
radius R calculate 
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The shift of the orders now becomes apparent: the in¬ 
tegral kernel contains the perturbation / explicitely in 
all terms. The calculation of the magnetic field up to 
order / 1 requires only the knowledge of the potential 
(and also G) up to order /°, i.e. for the spherical case. 
Thus up to first order the perturbation theories for the 
magnetic and electric case are disjoint. 

Please, note that (after partial integration) the integral 
vanishes if both / and V are axially symmetric (no 
$ dependence) in agreement to the general result that 
an axially symmetric source in an axially symmetric 
volume conductor induces no external magnetic field. 
The calculation of the total magnetic field from the 
radial component is straight forward: integrate along 
the radial direction to get the magnetic scalar poten¬ 
tial and then take the gradient resulting in an exactly 
curl-free solution. This step is trivial for an expan¬ 
sion in spherical harmonics (see next section) which 
are given in polar ccordinates. 


2.3 Algorithm 

To calculate the magnetic field up to first order we 
express V, G and / in series of spherical harmonics. 
The remaining integrals 
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2. For each sensor at position r, which measures B 
in direction % compute 
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Now, for the calculation of Bi, the magnetic field in 
the i.th sensor, we finally arrive at 
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where Vo denotes the gradient with respect to %. 

For notational convenience we have written the algo¬ 
rithm using complex numbers. However, in an ex¬ 
plicit computer implementation one can make use of 
the fact that V, G. and B are all real by splitting all 
terms into real and imaginary parts. 

The most important property of the theory is the 
sparseness of the ’matrices’ C nmpq ki. These fixed co¬ 
efficients assume non-vanishing values only if 


n = k + 2j with — ||>] #■$.< 2j < |p| - 1 
m = q ■ l . (14) 


Apart from the fact that sparseness saves computer 
time the first constraint ensures that for a finite 
parametrization of the surface the final series is an ex¬ 
pansion only in ro/r' where ro is the radius of the 
source and r 1 is the radius of the sensor (and not 



the volume conductor). Thus the series converges 
for sources arbitrily close to (or even on) the sur¬ 
face. In practice ro/r' is hardly larger than 0.6 even 
for very superficial sources. We found that up to 
ro/r' = 0.8 it is more than sufficient to sum the se¬ 
ries’ up to n, k < 20 ending up with a sparse and 
small (400 x 400) system matrix which has not to be 
inverted. 

3 Comparison for prolate spheroid 

We validated the theory for a spheroid for which an 
exact analytical solution is known. To calculate the 
solutions of the prolate spheroid we use our arbitrary- 
precision implementation [13] of the Cuffin/Cohen 
series expansion [3]. The relative differences of the 
solutions for 60 and 100 terms in the expansion were 
less (and usually far less) than 10 -4 for all sources 
used in this paper, which is sufficient to consider the 
solutions as ’exact’ for the comparisons to the differ¬ 
ent approximating solutions. 

The choice of the spherical approximation is not 
unique. We test two configuration: a) the best fit¬ 
ting sphere and b) the inner sphere (see Fig. 2). In 
Fig. 3 we show example fields for a deep dipole for 
which the spherical approximation completely fails. 
In Fig. 4 the relative deviation of the approximated 
solutions to the true solution is presented as a func¬ 
tion of dipole depth. In the lower panel we show 
the ratio of the errors arising from a measurement of 
the x-component and z-component of the magnetic 
field. While this ratio is essentially constant for ex¬ 
actly curl-free solutions (as is the presented perturba¬ 
tive solution) we eventually see a large dependence of 
the accuracy of BEM on the component under consid¬ 
eration. 

4 Discussion 

We solved the magnetic forward problem for realis¬ 
tic volume conductors analytically up to first order in 
the deviation of the realistic surface and a spherical 
approximation. This approach has several attractive 
properties 

1. The surface is approximated within perfectly 
smooth basis functions allowing to place a dipole 
arbitrarily close to the surface. 

2. For a finite surface parametrization the system 
matrix is sparse and its calculation requires no 



Figure 2: The prolate spheroid and two spherical ap¬ 
proximations. The approximation of the spheroid by 
spherical harmonics up to order P = 6 is practically 
exact and not included in the figure. The dipole loca¬ 
tion and orientation are varied in the analysis - the 
plot shows a typical example. Figure taken from [12] 



Figure 3: Calcutated fields for a dipole located at 
{xo,yo, zo) = (0, 0, —1) cm pointing in y-direction. 
Upper panel: exact solution (left), spherical approx¬ 
imation (middle), perturbative approximation (right). 
Lower panel: ideal correction to the sphere (left), cal¬ 
culated correction (middle), difference (ideal - calcu¬ 
lated) of the corrections (right). Figure taken from 
[ 12 ] 








Figure 4: Upper panel: Error of the forward cal¬ 
culation for various approximations of the forward 
calculation as a function of dipole height for the z- 
component of the magnetic field. Lower panel: Ratio 
of the corresponding error of the x-component and 
the z-component. The terms ’pert, a’ and ’pert, b’ 
refer to the perturbative expansions around a sphere 
of radius R = 10 cm and R = 9 cm, respectively. 
Figure taken from [12] 


inversion. 

3. Apart from unrealistic cases (like sources within 
the dewar) the system matrix can be chosen to be 
small (about 400 x 400). 

4. The solutions are exactly curl-free. 

However, certain requirements must be met to make 

the presented theory applicable in practice 

1. The surface can be reasonably approximated by 
a sphere (unlike a cigar). 

2. The surface is sufficiently smooth (unlike a golf 
ball). 
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